#!/bin/bash
set -e

log=.log
if [ ! -e "$log" ]; then
	:> $log
fi

while getopts  "s:p:" opts
do
        case  $opts  in
        s)
                        sample_name=$OPTARG
                        echo Sample name is: $sample_name 2>>$log 1>&2
                ;;
		p)
                        out_prefix=$OPTARG
                        echo out prefix is: $out_prefix 2>>$log 1>&2
                ;;
        esac
done
shift $(($OPTIND - 1))


if [ -z $sample_name ]; then
        sample_name=sap
        echo Sample name is: $sample_name 2>>$log 1>&2
fi


if [ $# -lt 2 ]; then
        echo "error... need args" 2>>$log 1>&2
        echo $0 [-s sample_name] [-p out_prefix] reads1.fq reads2.fq 2>>$log 1>&2
        exit 1
fi

if [ -z $out_prefix ]; then
	out_prefix=1
fi

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------
#-----------------------------------------------

###--------------- argument may be changed ---------------###


# genome_name=b37.fa
echo using $genome_name as the reference genome 2>>$log 1>&2
# genome_assembly=b37
echo genome assembly is: $genome_assembly 2>>$log 1>&2
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version 2>>$log 1>&2

# data_thread_num=8
# cpu_thread_num=4

# java_memory=16g
echo java memory: $java_memory 2>>$log 1>&2

# snpeff_db_version=GRCh37.75
echo snpeff database: $snpeff_db_version 2>>$log 1>&2

# header=\
# @HD\tVN:1.4\tGO:none\tSO:coordinate

read_group=\
@RG\\tID:${sample_name}\\tPL:ILLUMINA\\tSM:$sample_name
# For gatk, Each read group must contain the platform (PL) and sample (SM) tags.
# For the platform value, we currently support 454, LS454, Illumina, Solid, ABI_Solid, and CG (all case-insensitive).
# Each read in the file must be associated with exactly one read group.


# vcf_path=${data_path}/vcf/gatk/

# there is no space between "=" and "\"
# ref_genome=\
# ${data_path}/ref/b37/$genome_name

echo ref genome is: $ref_genome 2>>$log 1>&2

# reads_seq=\
# /mnt/lustre/users/medical001/Bre_rawdata/0914__hiseq_Qiagen_test/ccy/ccy_ACTTGA_L006.fq
# reads_seq_1=\
# /mnt/lustre-2/users/fei.yang/work/data/test/ccy.6.1.1.fq
# reads_seq_2=\
# /mnt/lustre-2/users/fei.yang/work/data/test/ccy.6.2.1.fq

# reads_seq_1=\
# /mnt/lustre/users/medical001/yangjm210/Sample_ZSL/ZSL_GATGAA_L003_R1_001.fq.gz
# reads_seq_2=\
# /mnt/lustre/users/medical001/yangjm210/Sample_ZSL/ZSL_GATGAA_L003_R2_001.fq.gz

reads_seq_1=$1
reads_seq_2=$2

echo reads are: $reads_seq_1 and $reads_seq_2 2>>$log 1>&2


###--------------- end argument may be changed ---------------###

###--------------- tools path ---------------###

# bwa=${tools_path}/bwa-0.7.12/bwa
# bwa=/share/apps/bwa-0.7.5a/bwa
# samtools=${tools_path}/samtools-1.2/samtools
# picard_path=${tools_path}/picard-tools-1.119/
# gatk_path=${tools_path}/GenomeAnalysisTK-1.6-9-g47df7bb/
# gatk3_path=${tools_path}
# gatk=${gatk3_path}/GenomeAnalysisTK.jar
# snpeff_path=${tools_path}/snpEff
# snpeff=${snpeff_path}/snpEff.jar
# snpsift=${snpeff_path}/SnpSift.jar
# fastuniq=${tools_path}/FastUniq/fastuniq
###--------------- end tools path ---------------###



# bwa index
# find $ref_genome.sa $ref_genome.amb $ref_genome.ann $ref_genome.pac $ref_genome.bwt
# if [ $? -ne 0 ]; then
	# echo Need to index ref genome ... 2>>$log 1>&2
	# $bwa index -a bwtsw $ref_genome
# else
	# echo Ref genome has been indexed 2>>$log 1>&2
# fi
# the outpt files are always the same dir with ref.fa
# output:
# hg19.fa.sa
# hg19.fa.amb
# hg19.fa.ann
# hg19.fa.pac
# hg19.fa.bwt

# for gatk
# find ${ref_genome}.fai
# if [ $? -ne 0 ]; then
	# $samtools faidx $ref_genome
# fi

# ref_genome_basename=`dirname $ref_genome`/`basename $ref_genome .fa`

# find ${ref_genome_basename}.dict
# if [ $? -ne 0 ]; then
	# java -Xmx$java_memory -jar ${picard_path}CreateSequenceDictionary.jar \
	# REFERENCE=$ref_genome \
	# OUTPUT=${ref_genome_basename}.dict \
	# GENOME_ASSEMBLY=$genome_assembly
# fi

 
echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo start bwa mem ... 2>>$log 1>&2
time \
$bwa mem -M -L 1000 -t 10 -R $read_group $ref_genome \
$reads_seq_1 \
$reads_seq_2 \
> $out_prefix.sam \
2>> $log


echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo mark primer 2>>$log 1>&2
${script_path}/mark_sam_qual0_start_end.sh $out_prefix.sam $out_prefix.mark.sam


echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo 'picard sam->bam sort index fixmateinformation' 2>>$log 1>&2
time \
java -Xmx$java_memory -jar ${picard_path}FixMateInformation.jar \
I=$out_prefix.mark.sam \
O=$out_prefix.sort.bam \
SORT_ORDER=coordinate \
CREATE_INDEX=true \
2>> $log

# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo 'picard sam->bam sort index fixmateinformation' 2>>$log 1>&2
# time \
# java -Xmx$java_memory -jar ${picard_path}FixMateInformation.jar \
# I=$out_prefix.sam \
# O=$out_prefix.sort.bam \
# SORT_ORDER=coordinate \
# CREATE_INDEX=true \
# 2>> $log




